# Convert dates to proper format
admissions$admittime = as.POSIXct(admissions$admittime)
admissions$dischtime = as.POSIXct(admissions$dischtime)

# Calculate length of stay
admissions$los_days = as.numeric(difftime(admissions$dischtime,
                                          admissions$admittime,
                                          units = 'days'))

# Remove invalid length of stay values
admissions = admissions %>% filter(los_days > 0 & los_days < 365)

cat('After filtering, admissions dataset contains', nrow(admissions), 'records\n')

# Identify readmissions within 30 days
admissions = admissions %>%
  arrange(subject_id, admittime) %>%
  group_by(subject_id) %>%
  mutate(
    next_admission = lead(admittime),
    days_to_readmit = as.numeric(difftime(next_admission, dischtime, units = 'days')),
    readmit_30day = ifelse(days_to_readmit <= 30 & !is.na(days_to_readmit), 1, 0)
  ) %>%
  ungroup()

# Filter to eligible discharges (exclude last admission per patient)
readmit_eligible = admissions %>%
  filter(!is.na(readmit_30day))

# Calculate overall readmission statistics
readmit_summary = readmit_eligible %>%
  summarise(
    total_discharges = n(),
    readmissions = sum(readmit_30day),
    readmission_rate = mean(readmit_30day) * 100,
    .groups = 'drop'
  )

# Display results
summary_table = data.frame(
  Metric = c('Total Eligible Discharges', '30-Day Readmissions', '30-day Readmission Rate (%)'),
  Value = c(readmit_summary$total_discharges,
            readmit_summary$readmissions,
            round(readmit_summary$readmission_rate, 2))
)

kable(summary_table, caption = '30-Day Readmission Summary Statistics')

cat('Eligible discharges for readmission analysis:', nrow(readmit_eligible), '\n')

# Merge with patient data for demographics
analysis_data = readmit_eligible %>%
  left_join(patients, by = 'subject_id') %>%
  mutate(age_at_adm = anchor_age)

# Age comparison by readmission status
age_stats = analysis_data %>%
  group_by(readmit_30day) %>%
  summarise(
    count = n(),
    mean_age = round(mean(age_at_adm, na.rm = TRUE), 1),
    median_age = round(median(age_at_adm, na.rm = TRUE), 1),
    sd_age = round(sd(age_at_adm, na.rm = TRUE), 1),
    .groups = 'drop'
  ) %>%
  mutate(readmission_status = ifelse(readmit_30day == 1, 'Readmitted', 'Not Readmitted'))

kable(age_stats[, c('readmission_status', 'count', 'mean_age', 'median_age', 'sd_age')],
      caption = 'Age Statistics by Readmission Status',
      col.names = c('Readmission Status', 'Count', 'Mean Age', 'Median Age', 'SD Age'))

# Length of stay comparison
los_stats <- analysis_data %>%
  group_by(readmit_30day) %>%
  summarise(
    count = n(),
    mean_los = round(mean(los_days, na.rm = TRUE), 1),
    median_los = round(median(los_days, na.rm = TRUE), 1),
    q75_los = round(quantile(los_days, 0.75, na.rm = TRUE), 1),
    max_los = round(max(los_days, na.rm = TRUE), 1),
    .groups = 'drop'
  ) %>%
  mutate(readmission_status = ifelse(readmit_30day == 1, "Readmitted", "Not Readmitted"))

kable(los_stats[, c("readmission_status", "count", "mean_los", "median_los", "q75_los")], 
      caption = "Length of Stay Statistics by Readmission Status",
      col.names = c("Readmission Status", "Count", "Mean LOS", "Median LOS", "75th Percentile LOS"))

# Basic overview of diagnoses data
cat("Diagnoses data contains", nrow(diagnoses), "records for", 
    length(unique(diagnoses$hadm_id)), "hospital admissions.\n")

cat("There are", length(unique(diagnoses$icd_code)), "unique ICD codes recorded.\n")

# Full integration with diagnosis dictionary
diagnoses_enhanced <- diagnoses %>%
  left_join(d_icd_diagnoses, by = c('icd_code' = 'icd_code', 'icd_version' = 'icd_version')) %>%
  mutate(
    # Use only long_title since short_title doesn't exist
    diagnosis_name = coalesce(long_title, paste("Unknown Diagnosis", icd_code)),
    # Create diagnosis categories based on ICD structure
    diagnosis_category = case_when(
      icd_version == 9 ~ case_when(
        substr(icd_code, 1, 3) %in% c("001", "002", "003") ~ "Infectious Diseases",
        substr(icd_code, 1, 1) == "1" ~ "Neoplasms",
        substr(icd_code, 1, 1) == "2" ~ "Endocrine/Metabolic",
        substr(icd_code, 1, 1) == "3" ~ "Blood Disorders",
        substr(icd_code, 1, 1) == "4" ~ "Circulatory System",
        substr(icd_code, 1, 1) == "5" ~ "Respiratory System",
        substr(icd_code, 1, 1) == "6" ~ "Digestive System",
        substr(icd_code, 1, 1) == "7" ~ "Genitourinary System",
        substr(icd_code, 1, 1) == "8" ~ "Injury/Poisoning",
        substr(icd_code, 1, 1) == "V" ~ "Supplementary Classification",
        TRUE ~ "Other ICD-9"
      ),
      icd_version == 10 ~ case_when(
        substr(icd_code, 1, 1) %in% c("A", "B") ~ "Infectious Diseases",
        substr(icd_code, 1, 1) == "C" ~ "Neoplasms",
        substr(icd_code, 1, 1) == "E" ~ "Endocrine/Metabolic",
        substr(icd_code, 1, 1) == "I" ~ "Circulatory System",
        substr(icd_code, 1, 1) == "J" ~ "Respiratory System",
        substr(icd_code, 1, 1) == "K" ~ "Digestive System",
        substr(icd_code, 1, 1) == "N" ~ "Genitourinary System",
        substr(icd_code, 1, 1) == "S" ~ "Injury/Poisoning",
        substr(icd_code, 1, 1) == "Z" ~ "Health Status Codes",
        TRUE ~ "Other ICD-10"
      ),
      TRUE ~ "Unknown Version"
    )
  )

# Check dictionary integration success
dict_integration_stats <- diagnoses_enhanced %>%
  summarise(
    total_records = n(),
    with_descriptions = sum(!is.na(long_title)),
    coverage_percent = round(sum(!is.na(long_title)) / n() * 100, 1),
    unique_categories = n_distinct(diagnosis_category)
  )

cat("\nDictionary Integration Results:\n")

cat("Total diagnosis records:", dict_integration_stats$total_records, "\n")

cat("Records with descriptions:", dict_integration_stats$with_descriptions, "\n")

cat("Dictionary coverage:", dict_integration_stats$coverage_percent, "%\n")

cat("Unique diagnosis categories:", dict_integration_stats$unique_categories, "\n")

# Analyze by diagnosis categories
category_summary <- diagnoses_enhanced %>%
  group_by(diagnosis_category) %>%
  summarise(
    total_cases = n(),
    unique_codes = n_distinct(icd_code),
    unique_patients = n_distinct(subject_id),
    unique_admissions = n_distinct(hadm_id),
    percent_of_total = round(n() / nrow(diagnoses_enhanced) * 100, 1),
    .groups = 'drop'
  ) %>%
  arrange(desc(total_cases))

kable(category_summary, caption = 'Diagnosis Distribution by Clinical Category',
      col.names = c('Clinical Category', 'Total Cases', 'Unique Codes', 'Unique Patients', 
                   'Unique Admissions', '% of Total'))

# Top diagnoses with full descriptions
icd_summary_enhanced <- diagnoses_enhanced %>%
  group_by(icd_code, icd_version, diagnosis_name, diagnosis_category) %>%
  summarise(
    frequency = n(),
    unique_patients = n_distinct(subject_id),
    unique_admissions = n_distinct(hadm_id),
    .groups = 'drop'
  ) %>%
  arrange(desc(frequency))

top_diagnoses_enhanced <- head(icd_summary_enhanced, 15) %>%
  mutate(
    # Truncate long descriptions for table readability
    diagnosis_display = ifelse(nchar(diagnosis_name) > 60, 
                              paste0(substr(diagnosis_name, 1, 57), "..."), 
                              diagnosis_name)
  )

kable(top_diagnoses_enhanced[, c('icd_code', 'icd_version', 'diagnosis_display', 'diagnosis_category', 'frequency', 'unique_patients')], 
      caption = 'Top 15 Most Frequent Diagnoses with Clinical Categories',
      col.names = c('ICD Code', 'Version', 'Diagnosis', 'Category', 'Frequency', 'Unique Patients'))

# Diagnoses per admission and readmission analysis using enhanced data
diagnoses_per_admission_enhanced = diagnoses_enhanced %>%
  group_by(hadm_id) %>%
  summarise(
    diagnoses_count = n(),
    unique_categories = n_distinct(diagnosis_category),
    has_circulatory = any(diagnosis_category == "Circulatory System"),
    has_respiratory = any(diagnosis_category == "Respiratory System"),
    has_endocrine = any(diagnosis_category == "Endocrine/Metabolic"),
    has_infectious = any(diagnosis_category == "Infectious Diseases"),
    .groups = 'drop'
  )

diagnoses_burden_table = diagnoses_per_admission_enhanced %>%
  summarise(
    mean_diagnoses = round(mean(diagnoses_count), 2),
    median_diagnoses = median(diagnoses_count),
    max_diagnoses = max(diagnoses_count),
    q75_diagnoses = quantile(diagnoses_count, 0.75),
    mean_categories = round(mean(unique_categories), 2)
  )

diagnoses_summary_table = data.frame(
  Metric = c('Mean Diagnoses per Admission', 'Median Diagnoses per Admission', 
             '75th Percentile Diagnoses', 'Max Diagnoses per Admission',
             'Mean Diagnostic Categories per Admission'),
  Value = c(diagnoses_burden_table$mean_diagnoses,
            diagnoses_burden_table$median_diagnoses,
            diagnoses_burden_table$q75_diagnoses,
            diagnoses_burden_table$max_diagnoses,
            diagnoses_burden_table$mean_categories)
)

kable(diagnoses_summary_table, caption = 'Enhanced Diagnoses Burden per Admission')

# Merge enhanced diagnoses data with analysis dataset
admission_diagnoses_enhanced = diagnoses_enhanced %>%
  group_by(hadm_id) %>%
  summarise(
    total_diagnoses = n(),
    unique_diagnosis_categories = n_distinct(diagnosis_category),
    primary_diagnosis_code = first(icd_code[seq_num == 1]),
    primary_diagnosis_name = first(diagnosis_name[seq_num == 1]),
    primary_diagnosis_category = first(diagnosis_category[seq_num == 1]),
    # High-risk diagnosis flags
    has_heart_failure = any(grepl("heart failure|cardiac failure", diagnosis_name, ignore.case = TRUE)),
    has_diabetes = any(grepl("diabetes", diagnosis_name, ignore.case = TRUE)),
    has_copd = any(grepl("chronic obstructive|copd|emphysema", diagnosis_name, ignore.case = TRUE)),
    has_pneumonia = any(grepl("pneumonia", diagnosis_name, ignore.case = TRUE)),
    .groups = 'drop'
  )

analysis_data = analysis_data %>%
  left_join(admission_diagnoses_enhanced, by = 'hadm_id')

# High-risk diagnoses and readmission analysis
high_risk_conditions <- analysis_data %>%
  summarise(
    heart_failure_readmit_rate = round(mean(readmit_30day[has_heart_failure == TRUE], na.rm = TRUE) * 100, 1),
    diabetes_readmit_rate = round(mean(readmit_30day[has_diabetes == TRUE], na.rm = TRUE) * 100, 1),
    copd_readmit_rate = round(mean(readmit_30day[has_copd == TRUE], na.rm = TRUE) * 100, 1),
    pneumonia_readmit_rate = round(mean(readmit_30day[has_pneumonia == TRUE], na.rm = TRUE) * 100, 1),
    overall_readmit_rate = round(mean(readmit_30day, na.rm = TRUE) * 100, 1)
  )

high_risk_table <- data.frame(
  Condition = c('Heart Failure', 'Diabetes', 'COPD', 'Pneumonia', 'Overall'),
  Readmission_Rate = c(high_risk_conditions$heart_failure_readmit_rate,
                       high_risk_conditions$diabetes_readmit_rate,
                       high_risk_conditions$copd_readmit_rate,
                       high_risk_conditions$pneumonia_readmit_rate,
                       high_risk_conditions$overall_readmit_rate)
)

kable(high_risk_table, caption = 'Readmission Rates by High-Risk Conditions',
      col.names = c('Condition', 'Readmission Rate (%)'))

# Primary diagnosis category analysis
primary_dx_category_analysis = analysis_data %>%
  group_by(primary_diagnosis_category) %>%
  summarise(
    total_cases = n(),
    readmissions = sum(readmit_30day, na.rm = TRUE),
    readmission_rate = round(mean(readmit_30day, na.rm = TRUE) * 100, 1),
    .groups = 'drop'
  ) %>%
  filter(total_cases >= 20 & !is.na(primary_diagnosis_category)) %>%
  arrange(desc(readmission_rate))

kable(primary_dx_category_analysis, caption = 'Readmission Rates by Primary Diagnosis Category',
      col.names = c('Primary Diagnosis Category', 'Total Cases', 'Readmissions', 'Readmission Rate (%)'))

# Enhanced diagnosis burden by readmission status
diagnosis_burden_enhanced = analysis_data %>%
  group_by(readmit_30day) %>%
  summarise(
    mean_diagnoses = round(mean(total_diagnoses, na.rm = TRUE), 2),
    median_diagnoses = median(total_diagnoses, na.rm = TRUE),
    mean_categories = round(mean(unique_diagnosis_categories, na.rm = TRUE), 2),
    q75_diagnoses = quantile(total_diagnoses, 0.75, na.rm = TRUE),
    .groups = 'drop'
  ) %>%
  mutate(readmission_status = ifelse(readmit_30day == 1, 'Readmitted', 'Not Readmitted'))

kable(diagnosis_burden_enhanced[, c('readmission_status', 'mean_diagnoses', 'median_diagnoses', 'mean_categories', 'q75_diagnoses')],
      caption = 'Enhanced Diagnosis Burden by Readmission Status',
      col.names = c('Readmission Status', 'Mean Diagnoses', 'Median Diagnoses', 'Mean Categories', '75th Percentile'))

cat('Prescriptions dataset contains', nrow(prescriptions), 'records\n')

cat('Unique patients:', length(unique(prescriptions$subject_id)), '\n')

cat('Unique admissions:', length(unique(prescriptions$hadm_id)), '\n')

cat('Unique medications:', length(unique(prescriptions$drug)), '\n')

# Enhanced prescription analysis with drug categorization
prescriptions_enhanced <- prescriptions %>%
  mutate(
    # Clean drug names for better categorization
    drug_clean = toupper(trimws(drug)),
    # Create therapeutic categories based on drug names
    drug_category = case_when(
        # IV Fluids and Electrolytes
        grepl("SODIUM CHLORIDE|DEXTROSE|LACTATED RINGERS|POTASSIUM CHLORIDE|MAGNESIUM SULFATE|CALCIUM GLUCONATE|STERILE WATER|GLUCOSE GEL|^NS$", drug_clean) ~ "IV Fluids/Electrolytes",
        
        # Pain Management
        grepl("HYDROMORPHONE|DILAUDID|MORPHINE|FENTANYL|OXYCODONE|TRAMADOL|ACETAMINOPHEN|IBUPROFEN|GABAPENTIN", drug_clean) ~ "Pain Management",
        
        # Anticoagulants (major category we missed)
        grepl("HEPARIN|WARFARIN|ASPIRIN", drug_clean) ~ "Anticoagulants",
        
        # CNS/Psychiatric Medications
        grepl("LORAZEPAM|ALPRAZOLAM|MIDAZOLAM|HALOPERIDOL", drug_clean) ~ "CNS Medications",
        
        # Proton Pump Inhibitors/GI
        grepl("PANTOPRAZOLE|OMEPRAZOLE|LANSOPRAZOLE|PROTONIX|SENNA|DOCUSATE|BISACODYL|POLYETHYLENE GLYCOL", drug_clean) ~ "GI Medications",
        
        # Cardiovascular (add statin)
        grepl("LISINOPRIL|METOPROLOL|ATENOLOL|AMLODIPINE|FUROSEMIDE|LASIX|ATORVASTATIN", drug_clean) ~ "Cardiovascular",
        
        # Corticosteroids
        grepl("PREDNISONE|PREDNISOLONE|METHYLPREDNISOLONE", drug_clean) ~ "Corticosteroids",
        
        # Anti-nausea
        grepl("ONDANSETRON|ZOFRAN", drug_clean) ~ "Anti-nausea",
        
        # Diabetes/Metabolic
        grepl("GLUCAGON|INSULIN|METFORMIN", drug_clean) ~ "Diabetes/Metabolic",
        
        # Antibiotics
        grepl("VANCOMYCIN|CIPROFLOXACIN|LEVOFLOXACIN|CEFTRIAXONE|AZITHROMYCIN", drug_clean) ~ "Antibiotics",
        
        # Administration/Container
        grepl("^BAG$|^VIAL$|^SW$|FLUSH", drug_clean) ~ "Administration/Container",
        
        TRUE ~ "Other/Unclassified"
      ),
    # High-risk medication flags
    high_risk_med = drug_category %in% c("Anticoagulants", "Opioid Analgesics", "CNS Medications"),
    # Route categorization
    route_category = case_when(
      grepl("IV|INTRAVENOUS", toupper(coalesce(route, ""))) ~ "Intravenous",
      grepl("PO|ORAL", toupper(coalesce(route, ""))) ~ "Oral",
      grepl("IM|INTRAMUSCULAR", toupper(coalesce(route, ""))) ~ "Intramuscular",
      grepl("SL|SUBLINGUAL", toupper(coalesce(route, ""))) ~ "Sublingual",
      grepl("TOP|TOPICAL", toupper(coalesce(route, ""))) ~ "Topical",
      TRUE ~ "Other/Unknown"
    )
  )

# Drug category analysis
drug_category_summary <- prescriptions_enhanced %>%
  group_by(drug_category) %>%
  summarise(
    total_prescriptions = n(),
    unique_drugs = n_distinct(drug),
    unique_patients = n_distinct(subject_id),
    unique_admissions = n_distinct(hadm_id),
    percent_of_total = round(n() / nrow(prescriptions_enhanced) * 100, 1),
    .groups = 'drop'
  ) %>%
  arrange(desc(total_prescriptions))

kable(drug_category_summary, caption = 'Prescription Distribution by Therapeutic Category',
      col.names = c('Therapeutic Category', 'Total Prescriptions', 'Unique Drugs', 
                   'Unique Patients', 'Unique Admissions', '% of Total'))

# Top medications with categories
medication_summary_enhanced = prescriptions_enhanced %>%
  group_by(drug, drug_category) %>%
  summarise(
    frequency = n(),
    unique_patients = n_distinct(subject_id),
    unique_admissions = n_distinct(hadm_id),
    .groups = 'drop'
  ) %>%
  arrange(desc(frequency)) %>%
  top_n(20, wt = frequency)

kable(medication_summary_enhanced, caption = 'Top 20 Medications by Frequency with Therapeutic Categories',
      col.names = c('Medication', 'Therapeutic Category', 'Total Prescriptions', 'Unique Patients', 'Unique Admissions'))

# Enhanced polypharmacy analysis with categories
medications_per_admission_enhanced = prescriptions_enhanced %>%
  group_by(hadm_id) %>%
  summarise(
    medication_count = n(),
    unique_medications = n_distinct(drug),
    unique_categories = n_distinct(drug_category),
    high_risk_med_count = sum(high_risk_med),
    has_diabetes_meds = any(drug_category == "Diabetes Medications"),
    has_cv_meds = any(drug_category == "Cardiovascular"),
    has_antibiotics = any(drug_category == "Antibiotics"),
    has_opioids = any(drug_category == "Opioid Analgesics"),
    .groups = 'drop'
  )

polypharmacy_enhanced_stats = medications_per_admission_enhanced %>%
  summarise(
    mean_medications = round(mean(medication_count), 1),
    median_medications = median(medication_count),
    max_medications = max(medication_count),
    q75_medications = quantile(medication_count, 0.75),
    mean_categories = round(mean(unique_categories), 1),
    median_categories = median(unique_categories),
    mean_high_risk = round(mean(high_risk_med_count), 1)
  )

polypharmacy_enhanced_table = data.frame(
  Metric = c('Mean Medications per Admission', 'Median Medications per Admission', 
             '75th Percentile Medications', 'Max Medications per Admission',
             'Mean Therapeutic Categories', 'Median Therapeutic Categories',
             'Mean High-Risk Medications'),
  Value = c(polypharmacy_enhanced_stats$mean_medications,
            polypharmacy_enhanced_stats$median_medications,
            polypharmacy_enhanced_stats$q75_medications,
            polypharmacy_enhanced_stats$max_medications,
            polypharmacy_enhanced_stats$mean_categories,
            polypharmacy_enhanced_stats$median_categories,
            polypharmacy_enhanced_stats$mean_high_risk)
)

kable(polypharmacy_enhanced_table, caption = 'Enhanced Polypharmacy Metrics per Admission')

# Remove existing medication variables before joining
analysis_data <- analysis_data %>%
  select(-any_of(c("medication_count", "unique_medications", "unique_categories", 
                   "high_risk_med_count", "has_diabetes_meds", "has_cv_meds", 
                   "has_antibiotics", "has_opioids")))

# Merge enhanced prescription data with analysis dataset AND clean in one step
analysis_data = analysis_data %>%
  left_join(medications_per_admission_enhanced, by = 'hadm_id') %>%
  mutate(
    medication_count = ifelse(is.na(medication_count), 0, medication_count),
    unique_medications = ifelse(is.na(unique_medications), 0, unique_medications),
    unique_categories = ifelse(is.na(unique_categories), 0, unique_categories),
    high_risk_med_count = ifelse(is.na(high_risk_med_count), 0, high_risk_med_count),
    has_diabetes_meds = ifelse(is.na(has_diabetes_meds), FALSE, has_diabetes_meds),
    has_cv_meds = ifelse(is.na(has_cv_meds), FALSE, has_cv_meds),
    has_antibiotics = ifelse(is.na(has_antibiotics), FALSE, has_antibiotics),
    has_opioids = ifelse(is.na(has_opioids), FALSE, has_opioids)
  )

# High-risk medication analysis and readmission
high_risk_med_analysis <- analysis_data %>%
  summarise(
    diabetes_meds_readmit_rate = round(mean(readmit_30day[has_diabetes_meds == TRUE], na.rm = TRUE) * 100, 1),
    cv_meds_readmit_rate = round(mean(readmit_30day[has_cv_meds == TRUE], na.rm = TRUE) * 100, 1),
    antibiotics_readmit_rate = round(mean(readmit_30day[has_antibiotics == TRUE], na.rm = TRUE) * 100, 1),
    opioids_readmit_rate = round(mean(readmit_30day[has_opioids == TRUE], na.rm = TRUE) * 100, 1),
    overall_readmit_rate = round(mean(readmit_30day, na.rm = TRUE) * 100, 1)
  )

med_risk_table <- data.frame(
  Medication_Category = c('Diabetes Medications', 'Cardiovascular', 'Antibiotics', 'Opioids', 'Overall'),
  Readmission_Rate = c(high_risk_med_analysis$diabetes_meds_readmit_rate,
                       high_risk_med_analysis$cv_meds_readmit_rate,
                       high_risk_med_analysis$antibiotics_readmit_rate,
                       high_risk_med_analysis$opioids_readmit_rate,
                       high_risk_med_analysis$overall_readmit_rate)
)

kable(med_risk_table, caption = 'Readmission Rates by Medication Categories',
      col.names = c('Medication Category', 'Readmission Rate (%)'))

# Enhanced prescription burden by readmission status
prescription_readmission_enhanced = analysis_data %>%
  group_by(readmit_30day) %>%
  summarise(
    mean_prescriptions = round(mean(medication_count, na.rm = TRUE), 2),
    median_prescriptions = median(medication_count, na.rm = TRUE),
    mean_categories = round(mean(unique_categories, na.rm = TRUE), 2),
    mean_high_risk = round(mean(high_risk_med_count, na.rm = TRUE), 2),
    q75_prescriptions = quantile(medication_count, 0.75, na.rm = TRUE),
    .groups = 'drop'
  ) %>%
  mutate(readmission_status = ifelse(readmit_30day == 1, 'Readmitted', 'Not Readmitted'))

kable(prescription_readmission_enhanced[, c('readmission_status', 'mean_prescriptions', 'median_prescriptions', 'mean_categories', 'mean_high_risk')],
      caption = 'Enhanced Prescription Burden by Readmission Status',
      col.names = c('Readmission Status', 'Mean Prescriptions', 'Median Prescriptions', 'Mean Categories', 'Mean High-Risk Meds'))

cat('\nEnhanced prescription analysis complete with full integration into analysis_data.\n')

cat('Lab events dataset contains', nrow(labevents), 'records\n')

cat('Unique patients:', length(unique(labevents$subject_id)), '\n')

cat('Unique admissions:', length(unique(labevents$hadm_id)), '\n')

cat('Unique lab tests:', length(unique(labevents$itemid)), '\n')

# Memory-efficient approach: Skip individual lab enhancement, go straight to aggregates
labs_per_admission_enhanced = labevents %>%
  left_join(d_labitems, by = 'itemid') %>%
  group_by(hadm_id) %>%
  summarise(
    total_lab_tests = n(),
    unique_lab_types = n_distinct(itemid),
    unique_categories = n_distinct(category, na.rm = TRUE),
    # Count specific important lab categories
    chemistry_tests = sum(grepl("chemistry", category, ignore.case = TRUE), na.rm = TRUE),
    hematology_tests = sum(grepl("hematology|blood", category, ignore.case = TRUE), na.rm = TRUE),
    has_labs = TRUE,
    .groups = 'drop'
  )

# Basic lab statistics
lab_burden_stats = labs_per_admission_enhanced %>%
  summarise(
    mean_tests = round(mean(total_lab_tests), 1),
    median_tests = median(total_lab_tests),
    max_tests = max(total_lab_tests),
    q75_tests = quantile(total_lab_tests, 0.75),
    mean_categories = round(mean(unique_categories), 1),
    median_categories = median(unique_categories)
  )

lab_burden_table = data.frame(
  Metric = c('Mean Lab Tests per Admission', 'Median Lab Tests per Admission', 
             '75th Percentile Lab Tests', 'Max Lab Tests per Admission',
             'Mean Lab Categories per Admission', 'Median Lab Categories per Admission'),
  Value = c(lab_burden_stats$mean_tests,
            lab_burden_stats$median_tests,
            lab_burden_stats$q75_tests,
            lab_burden_stats$max_tests,
            lab_burden_stats$mean_categories,
            lab_burden_stats$median_categories)
)

kable(lab_burden_table, caption = 'Laboratory Testing Burden per Admission')

# Remove existing lab variables before joining (for reproducibility)
analysis_data <- analysis_data %>%
  select(-any_of(c("total_lab_tests", "unique_lab_types", "unique_categories", 
                   "chemistry_tests", "hematology_tests", "had_labs", "has_labs")))

# Merge with analysis dataset
analysis_data = analysis_data %>%
  left_join(labs_per_admission_enhanced, by = 'hadm_id') %>%
  mutate(
    total_lab_tests = ifelse(is.na(total_lab_tests), 0, total_lab_tests),
    unique_lab_types = ifelse(is.na(unique_lab_types), 0, unique_lab_types),
    unique_categories = ifelse(is.na(unique_categories), 0, unique_categories),
    chemistry_tests = ifelse(is.na(chemistry_tests), 0, chemistry_tests),
    hematology_tests = ifelse(is.na(hematology_tests), 0, hematology_tests),
    has_labs = ifelse(is.na(has_labs), FALSE, has_labs)
  )

# Lab burden by readmission status
lab_readmission_analysis = analysis_data %>%
  group_by(readmit_30day) %>%
  summarise(
    mean_lab_tests = round(mean(total_lab_tests, na.rm = TRUE), 1),
    median_lab_tests = median(total_lab_tests, na.rm = TRUE),
    mean_categories = round(mean(unique_categories, na.rm = TRUE), 1),
    percent_with_labs = round(mean(has_labs, na.rm = TRUE) * 100, 1),
    q75_lab_tests = quantile(total_lab_tests, 0.75, na.rm = TRUE),
    .groups = 'drop'
  ) %>%
  mutate(readmission_status = ifelse(readmit_30day == 1, 'Readmitted', 'Not Readmitted'))

kable(lab_readmission_analysis[, c('readmission_status', 'mean_lab_tests', 'median_lab_tests', 
                                   'mean_categories', 'percent_with_labs')],
      caption = 'Laboratory Testing Burden by Readmission Status',
      col.names = c('Readmission Status', 'Mean Lab Tests', 'Median Lab Tests', 
                   'Mean Categories', '% with Labs'))

cat('\nLab data coverage:\n')

cat('Admissions with lab data:', sum(analysis_data$total_lab_tests > 0), '/', nrow(analysis_data), '\n')

cat('Coverage:', round(sum(analysis_data$total_lab_tests > 0) / nrow(analysis_data) * 100, 1), '%\n')

# Simplified Lab Visualization: Lab testing volume per admission
ggplot(labs_per_admission_enhanced, aes(x = total_lab_tests)) +
  geom_histogram(bins = 50, fill = 'orange', alpha = 0.7) +
  labs(title = 'Distribution of Total Lab Tests per Admission',
       x = 'Number of Lab Tests',
       y = 'Number of Admissions') +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

# Lab categories distribution
ggplot(labs_per_admission_enhanced, aes(x = unique_categories)) +
  geom_histogram(bins = 20, fill = 'darkgreen', alpha = 0.7) +
  labs(title = 'Distribution of Lab Categories per Admission',
       x = 'Number of Lab Categories',
       y = 'Number of Admissions') +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

# Remove existing lab variables before joining (for reproducibility)
analysis_data <- analysis_data %>%
  select(-any_of(c("total_lab_tests", "unique_lab_types", "unique_categories", 
                   "chemistry_tests", "hematology_tests", "had_labs", "has_labs")))

# Merge lab data with analysis dataset
analysis_data = analysis_data %>%
  left_join(labs_per_admission_enhanced, by = 'hadm_id') %>%
  mutate(
    total_lab_tests = ifelse(is.na(total_lab_tests), 0, total_lab_tests),
    unique_lab_types = ifelse(is.na(unique_lab_types), 0, unique_lab_types),
    unique_categories = ifelse(is.na(unique_categories), 0, unique_categories),
    chemistry_tests = ifelse(is.na(chemistry_tests), 0, chemistry_tests),
    hematology_tests = ifelse(is.na(hematology_tests), 0, hematology_tests),
    has_labs = ifelse(is.na(has_labs), FALSE, has_labs)
  )

# Lab burden by readmission status
lab_readmission_analysis = analysis_data %>%
  group_by(readmit_30day) %>%
  summarise(
    mean_lab_tests = round(mean(total_lab_tests, na.rm = TRUE), 1),
    median_lab_tests = median(total_lab_tests, na.rm = TRUE),
    q75_lab_tests = quantile(total_lab_tests, 0.75, na.rm = TRUE),
    percent_with_labs = round(mean(has_labs, na.rm = TRUE) * 100, 1),
    mean_unique_types = round(mean(unique_lab_types, na.rm = TRUE), 1),
    .groups = 'drop'
  ) %>%
  mutate(readmission_status = ifelse(readmit_30day == 1, 'Readmitted', 'Not Readmitted'))

kable(lab_readmission_analysis[, c('readmission_status', 'mean_lab_tests', 'median_lab_tests', 
                                   'q75_lab_tests', 'percent_with_labs', 'mean_unique_types')],
      caption = 'Laboratory Testing Burden by Readmission Status',
      col.names = c('Readmission Status', 'Mean Lab Tests', 'Median Lab Tests', 
                   '75th Percentile', '% with Labs', 'Mean Unique Types'))

# Statistical test for difference in lab burden
lab_test_result <- t.test(total_lab_tests ~ readmit_30day, data = analysis_data)
cat('\nT-test for difference in lab test burden between groups:\n')

cat('t =', round(lab_test_result$statistic, 3), ', p-value =', 
    format(lab_test_result$p.value, scientific = TRUE, digits = 3), '\n')

# Calculate coverage
cat('\nLab data coverage:\n')

cat('Admissions with lab data:', sum(analysis_data$total_lab_tests > 0), '/', nrow(analysis_data), '\n')

cat('Coverage:', round(sum(analysis_data$total_lab_tests > 0) / nrow(analysis_data) * 100, 1), '%\n')

# Skip individual lab analysis to avoid memory issues
cat('\n=== INDIVIDUAL LAB ANALYSIS SKIPPED ===\n')

cat('(Individual lab analysis removed to prevent memory limit errors)\n')

cat('Procedures dataset contains', nrow(procedures), 'records\n')

cat('Unique patients:', length(unique(procedures$subject_id)), '\n')

cat('Unique admissions:', length(unique(procedures$hadm_id)), '\n')

cat('Unique procedures:', length(unique(procedures$icd_code)), '\n')

# Enhanced procedures analysis with full dictionary integration
procedures_enhanced <- procedures %>%
  left_join(d_icd_procedures, by = c('icd_code' = 'icd_code', 'icd_version' = 'icd_version')) %>%
  mutate(
    # Create meaningful procedure name
    procedure_name = coalesce(long_title, paste("Unknown Procedure", icd_code)),
    
    # Create procedure categories based on ICD structure and clinical knowledge
    procedure_category = case_when(
      icd_version == 9 ~ case_when(
        substr(icd_code, 1, 2) %in% c("00", "01") ~ "CNS Procedures",
        substr(icd_code, 1, 2) %in% c("02", "03", "04", "05") ~ "Endocrine/Eye/ENT Procedures",
        substr(icd_code, 1, 2) %in% c("06", "07", "08", "09") ~ "Respiratory Procedures",
        substr(icd_code, 1, 2) %in% c("35", "36", "37", "38", "39") ~ "Cardiovascular Procedures",
        substr(icd_code, 1, 2) %in% c("42", "43", "44", "45", "46") ~ "Digestive System Procedures",
        substr(icd_code, 1, 2) %in% c("50", "51", "52", "53", "54") ~ "Urogenital Procedures",
        substr(icd_code, 1, 2) %in% c("77", "78", "79", "80", "81") ~ "Musculoskeletal Procedures",
        substr(icd_code, 1, 2) %in% c("86", "87", "88", "89") ~ "Integumentary/Diagnostic Procedures",
        substr(icd_code, 1, 2) %in% c("92", "93", "94", "95") ~ "Therapeutic Procedures",
        substr(icd_code, 1, 2) %in% c("96", "97", "98", "99") ~ "Miscellaneous Procedures",
        TRUE ~ "Other ICD-9 Procedures"
      ),
      icd_version == 10 ~ case_when(
        substr(icd_code, 1, 1) == "0" ~ case_when(
          substr(icd_code, 2, 2) %in% c("0", "1", "2") ~ "CNS Procedures",
          substr(icd_code, 2, 2) %in% c("3", "4", "5") ~ "Cardiovascular Procedures", 
          substr(icd_code, 2, 2) %in% c("6", "7", "8") ~ "Respiratory/Digestive Procedures",
          substr(icd_code, 2, 2) == "9" ~ "Reproductive Procedures",
          TRUE ~ "Other System Procedures"
        ),
        substr(icd_code, 1, 1) == "3" ~ "Administration Procedures",
        substr(icd_code, 1, 1) == "4" ~ "Measurement/Monitoring",
        substr(icd_code, 1, 1) == "5" ~ "Extracorporeal Assistance",
        substr(icd_code, 1, 1) == "6" ~ "Extracorporeal Therapies",
        substr(icd_code, 1, 1) == "7" ~ "Osteopathic Procedures",
        substr(icd_code, 1, 1) == "8" ~ "Other Procedures",
        substr(icd_code, 1, 1) == "9" ~ "Chiropractic Procedures",
        substr(icd_code, 1, 1) == "B" ~ "Imaging Procedures",
        substr(icd_code, 1, 1) == "C" ~ "Nuclear Medicine",
        substr(icd_code, 1, 1) == "D" ~ "Radiation Therapy",
        substr(icd_code, 1, 1) == "F" ~ "Physical Rehabilitation",
        substr(icd_code, 1, 1) == "G" ~ "Mental Health Procedures",
          TRUE ~ "Other ICD-10 Procedures"
      ),
      TRUE ~ "Unknown Version"
    ),
    # High-risk procedure flags based on clinical knowledge
    high_risk_procedure = procedure_category %in% c("CNS Procedures", "Cardiovascular Procedures") |
                         grepl("surgery|transplant|bypass|catheter", procedure_name, ignore.case = TRUE),
    # Invasive vs non-invasive classification
    invasive_procedure = !procedure_category %in% c("Imaging Procedures", "Nuclear Medicine", 
                                                   "Measurement/Monitoring", "Physical Rehabilitation") &
                        !grepl("consultation|examination|evaluation", procedure_name, ignore.case = TRUE)
  )

# Check dictionary integration success for procedures
proc_dict_integration_stats <- procedures_enhanced %>%
  summarise(
    total_records = n(),
    with_descriptions = sum(!is.na(long_title)),
    coverage_percent = round(sum(!is.na(long_title)) / n() * 100, 1),
    unique_categories = n_distinct(procedure_category)
  )

cat("\nProcedure Dictionary Integration Results:\n")

cat("Total procedure records:", proc_dict_integration_stats$total_records, "\n")

cat("Records with descriptions:", proc_dict_integration_stats$with_descriptions, "\n")

cat("Dictionary coverage:", proc_dict_integration_stats$coverage_percent, "%\n")

cat("Unique procedure categories:", proc_dict_integration_stats$unique_categories, "\n")

# Procedure category analysis
proc_category_summary <- procedures_enhanced %>%
  group_by(procedure_category) %>%
  summarise(
    total_procedures = n(),
    unique_codes = n_distinct(icd_code),
    unique_patients = n_distinct(subject_id),
    unique_admissions = n_distinct(hadm_id),
    percent_of_total = round(n() / nrow(procedures_enhanced) * 100, 1),
    .groups = 'drop'
  ) %>%
  arrange(desc(total_procedures))

kable(proc_category_summary, caption = 'Procedure Distribution by Clinical Category',
      col.names = c('Clinical Category', 'Total Procedures', 'Unique Codes', 'Unique Patients', 
                   'Unique Admissions', '% of Total'))

# Top procedures with full descriptions
procedure_summary_enhanced <- procedures_enhanced %>%
  group_by(icd_code, icd_version, procedure_name, procedure_category, high_risk_procedure) %>%
  summarise(
    frequency = n(),
    unique_patients = n_distinct(subject_id),
    unique_admissions = n_distinct(hadm_id),
    .groups = 'drop'
  ) %>%
  arrange(desc(frequency))

top_procedures_enhanced <- head(procedure_summary_enhanced, 15) %>%
  mutate(
    # Truncate long descriptions for table readability
    procedure_display = ifelse(nchar(procedure_name) > 50, 
                              paste0(substr(procedure_name, 1, 47), "..."), 
                              procedure_name),
    risk_level = ifelse(high_risk_procedure, "High Risk", "Standard")
  )

kable(top_procedures_enhanced[, c('icd_code', 'icd_version', 'procedure_display', 'procedure_category', 'risk_level', 'frequency')], 
      caption = 'Top 15 Most Frequent Procedures with Clinical Categories',
      col.names = c('ICD Code', 'Version', 'Procedure', 'Category', 'Risk Level', 'Frequency'))

# Enhanced procedures per admission analysis
procedures_per_admission_enhanced = procedures_enhanced %>%
  group_by(hadm_id) %>%
  summarise(
    procedure_count = n(),
    unique_procedures = n_distinct(icd_code),
    unique_categories = n_distinct(procedure_category),
    high_risk_procedure_count = sum(high_risk_procedure),
    invasive_procedure_count = sum(invasive_procedure),
    has_cardiovascular_proc = any(procedure_category == "Cardiovascular Procedures"),
    has_respiratory_proc = any(procedure_category == "Respiratory Procedures" | procedure_category == "Respiratory/Digestive Procedures"),
    has_cns_proc = any(procedure_category == "CNS Procedures"),
    has_imaging = any(procedure_category == "Imaging Procedures"),
    .groups = 'drop'
  )

procedure_burden_enhanced_stats = procedures_per_admission_enhanced %>%
  summarise(
    mean_procedures = round(mean(procedure_count), 1),
    median_procedures = round(median(procedure_count), 1),
    max_procedures = max(procedure_count),
    q75_procedures = round(quantile(procedure_count, 0.75), 1),
    mean_categories = round(mean(unique_categories), 1),
    mean_high_risk = round(mean(high_risk_procedure_count), 1),
    mean_invasive = round(mean(invasive_procedure_count), 1)
  )

procedure_burden_enhanced_table <- data.frame(
  Metric = c('Mean Procedures per Admission', 'Median Procedures per Admission', 
             'Max Procedures per Admission', '75th Percentile Procedures',
             'Mean Procedure Categories', 'Mean High-Risk Procedures', 'Mean Invasive Procedures'),
  Value = c(procedure_burden_enhanced_stats$mean_procedures, 
            procedure_burden_enhanced_stats$median_procedures,
            procedure_burden_enhanced_stats$max_procedures,
            procedure_burden_enhanced_stats$q75_procedures,
            procedure_burden_enhanced_stats$mean_categories,
            procedure_burden_enhanced_stats$mean_high_risk,
            procedure_burden_enhanced_stats$mean_invasive)
)

kable(procedure_burden_enhanced_table, caption = 'Enhanced Procedure Burden per Admission')

# Remove existing procedure variables before joining (for reproducibility)
analysis_data <- analysis_data %>%
  select(-any_of(c("procedure_count", "unique_procedures", "unique_categories", 
                   "high_risk_procedure_count", "invasive_procedure_count", 
                   "has_cardiovascular_proc", "has_respiratory_proc", "has_cns_proc", "has_imaging")))

# Merge enhanced procedure data with analysis dataset
analysis_data = analysis_data %>%
  left_join(procedures_per_admission_enhanced, by = 'hadm_id') %>%
  mutate(
    procedure_count = ifelse(is.na(procedure_count), 0, procedure_count),
    unique_procedures = ifelse(is.na(unique_procedures), 0, unique_procedures),
    unique_categories = ifelse(is.na(unique_categories), 0, unique_categories),
    high_risk_procedure_count = ifelse(is.na(high_risk_procedure_count), 0, high_risk_procedure_count),
    invasive_procedure_count = ifelse(is.na(invasive_procedure_count), 0, invasive_procedure_count),
    has_cardiovascular_proc = ifelse(is.na(has_cardiovascular_proc), FALSE, has_cardiovascular_proc),
    has_respiratory_proc = ifelse(is.na(has_respiratory_proc), FALSE, has_respiratory_proc),
    has_cns_proc = ifelse(is.na(has_cns_proc), FALSE, has_cns_proc),
    has_imaging = ifelse(is.na(has_imaging), FALSE, has_imaging)
  )

# High-risk procedure analysis and readmission
high_risk_proc_analysis <- analysis_data %>%
  summarise(
    cv_proc_readmit_rate = round(mean(readmit_30day[has_cardiovascular_proc == TRUE], na.rm = TRUE) * 100, 1),
    respiratory_proc_readmit_rate = round(mean(readmit_30day[has_respiratory_proc == TRUE], na.rm = TRUE) * 100, 1),
    cns_proc_readmit_rate = round(mean(readmit_30day[has_cns_proc == TRUE], na.rm = TRUE) * 100, 1),
    imaging_readmit_rate = round(mean(readmit_30day[has_imaging == TRUE], na.rm = TRUE) * 100, 1),
    overall_readmit_rate = round(mean(readmit_30day, na.rm = TRUE) * 100, 1)
  )

proc_risk_table <- data.frame(
  Procedure_Category = c('Cardiovascular', 'Respiratory', 'CNS', 'Imaging', 'Overall'),
  Readmission_Rate = c(high_risk_proc_analysis$cv_proc_readmit_rate,
                       high_risk_proc_analysis$respiratory_proc_readmit_rate,
                       high_risk_proc_analysis$cns_proc_readmit_rate,
                       high_risk_proc_analysis$imaging_readmit_rate,
                       high_risk_proc_analysis$overall_readmit_rate)
)

kable(proc_risk_table, caption = 'Readmission Rates by Procedure Categories',
      col.names = c('Procedure Category', 'Readmission Rate (%)'))

# Enhanced procedure burden by readmission status
procedure_readmission_enhanced = analysis_data %>%
  group_by(readmit_30day) %>%
  summarise(
    mean_procedures = round(mean(procedure_count, na.rm = TRUE), 2),
    median_procedures = median(procedure_count, na.rm = TRUE),
    mean_categories = round(mean(unique_categories, na.rm = TRUE), 2),
    mean_high_risk = round(mean(high_risk_procedure_count, na.rm = TRUE), 2),
    mean_invasive = round(mean(invasive_procedure_count, na.rm = TRUE), 2),
    q75_procedures = quantile(procedure_count, 0.75, na.rm = TRUE),
    .groups = 'drop'
  ) %>%
  mutate(readmission_status = ifelse(readmit_30day == 1, 'Readmitted', 'Not Readmitted'))

kable(procedure_readmission_enhanced[, c('readmission_status', 'mean_procedures', 'median_procedures', 'mean_categories', 'mean_high_risk', 'mean_invasive')],
      caption = 'Enhanced Procedure Burden by Readmission Status',
      col.names = c('Readmission Status', 'Mean Procedures', 'Median Procedures', 'Mean Categories', 'Mean High-Risk', 'Mean Invasive'))

# Readmission rate by admission type
admission_type_summary = analysis_data %>%
  group_by(admission_type) %>%
  summarise(
    total = n(),
    readmissions = sum(readmit_30day),
    rate = (readmissions / total) * 100,
    .groups = 'drop'
  ) %>%
  filter(total >= 10) # Only show types with sufficient sample size

ggplot(admission_type_summary, aes(x = reorder(admission_type, rate), y = rate)) +
  geom_col(fill = 'steelblue', alpha = 0.7) + 
  geom_text(aes(label = paste0(round(rate, 1), '%\n(', readmissions, '/', total, ')')),
            hjust = -0.1, size = 3) +
  coord_flip() +
  labs(title = '30-Day Readmission Rate by Admission Type',
       x = 'Admission Type',
       y = 'Readmission Rate (%)') +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

# Age distribution by readmission status
ggplot(analysis_data, aes(x = age_at_adm, fill = factor(readmit_30day))) +
  geom_histogram(position = 'identity', alpha = 0.6, bins = 30) +
  scale_fill_manual(values = c('0' = 'lightblue', '1' = 'salmon'),
                    labels = c('0' = 'Not Readmitted', '1' = 'Readmitted'),
                    name = 'Readmission Status') +
  labs(title = 'Age Distribution by Readmission Status',
       x = 'Age at Admission',
       y = 'Count') +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

# Clinical burden comparison
clinical_burden_plot_data <- analysis_data %>%
  select(readmit_30day, total_diagnoses, medication_count, procedure_count, total_lab_tests) %>%
  pivot_longer(cols = c(total_diagnoses, medication_count, procedure_count, total_lab_tests),
               names_to = 'clinical_measure',
               values_to = 'count') %>%
  mutate(
    clinical_measure = case_when(
      clinical_measure == 'total_diagnoses' ~ 'Diagnoses',
      clinical_measure == 'medication_count' ~ 'Prescriptions',
      clinical_measure == 'procedure_count' ~ 'Procedures',
      clinical_measure == 'total_lab_tests' ~ 'Lab Tests'
    ),
    readmission_status = ifelse(readmit_30day == 1, 'Readmitted', 'Not Readmitted')
  )

ggplot(clinical_burden_plot_data, aes(x = factor(readmit_30day), y = count, fill = factor(readmit_30day))) +
  geom_boxplot(alpha = 0.7) +
  facet_wrap(~ clinical_measure, scales = 'free_y') +
  scale_fill_manual(values = c('0' = 'lightblue', '1' = 'salmon'),
                    labels = c('0' = 'Not Readmitted', '1' = 'Readmitted'),
                    name = 'Readmission Status') +
  labs(title = 'Clinical Burden by Readmission Status',
       x = 'Readmission Status',
       y = 'Count') +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

cat("DEBUGGING:", ncol(analysis_data), "\n")

# Assess missing data patterns for key variables
key_vars <- c('subject_id', 'hadm_id', 'readmit_30day', 'age_at_adm', 'gender', 
              'race', 'insurance', 'admission_type', 'los_days', 'total_diagnoses',
              'unique_diagnosis_categories', 'medication_count', 'unique_categories', 
              'procedure_count', 'total_lab_tests', 'unique_lab_types')

completeness_summary <- analysis_data %>%
  select(all_of(key_vars)) %>%
  summarise(across(everything(), ~ sum(!is.na(.)))) %>%
  pivot_longer(everything(), names_to = 'variable', values_to = 'complete_count') %>%
  mutate(
    total_cases = nrow(analysis_data),
    completeness_percent = round(complete_count / total_cases * 100, 1),
    missing_count = total_cases - complete_count,
    missing_percent = round(missing_count / total_cases * 100, 1)
  ) %>%
  arrange(desc(completeness_percent))

kable(completeness_summary[, c('variable', 'complete_count', 'missing_count', 'completeness_percent', 'missing_percent')],
      caption = 'Data Completeness Summary for Key Variables',
      col.names = c('Variable', 'Complete Count', 'Missing Count', 'Completeness %', 'Missing %'))

# Calculate final analysis dataset characteristics
cat('\n=== FINAL ANALYSIS DATASET SUMMARY ===\n')

cat('Total admissions in analysis dataset:', nrow(analysis_data), '\n')

cat('Patients with readmissions:', sum(analysis_data$readmit_30day, na.rm = TRUE), '\n')

cat('Overall readmission rate:', round(mean(analysis_data$readmit_30day, na.rm = TRUE) * 100, 2), '%\n')

cat('Coverage rates:\n')

cat('- Diagnoses data:', round(sum(analysis_data$total_diagnoses > 0, na.rm = TRUE) / nrow(analysis_data) * 100, 1), '%\n')

cat('- Prescriptions data:', round(sum(analysis_data$medication_count > 0, na.rm = TRUE) / nrow(analysis_data) * 100, 1), '%\n')

cat('- Procedures data:', round(sum(analysis_data$procedure_count > 0, na.rm = TRUE) / nrow(analysis_data) * 100, 1), '%\n')

cat('- Lab data:', round(sum(analysis_data$total_lab_tests > 0, na.rm = TRUE) / nrow(analysis_data) * 100, 1), '%\n')

cat('- High-risk conditions identified:', round(sum(analysis_data$has_heart_failure | analysis_data$has_diabetes | analysis_data$has_copd, na.rm = TRUE) / nrow(analysis_data) * 100, 1), '%\n')

# Charlson Comorbidity Index calculation -- Healthcare standard method
# Based on ICD-9 and ICD-10 codes

# Charlson Comorbidity Index with proper NA handling
charlson_components <- analysis_data %>%
  mutate(
    # Cardiovascular conditions
    charlson_mi = ifelse(!is.na(primary_diagnosis_name) & grepl("myocardial infarction|acute mi", primary_diagnosis_name, ignore.case = TRUE), 1, 0),
    charlson_chf = ifelse(has_heart_failure, 1, 0),
    charlson_pvd = ifelse(!is.na(primary_diagnosis_name) & grepl("peripheral vascular|arterial disease", primary_diagnosis_name, ignore.case = TRUE), 1, 0),
    charlson_cvd = ifelse(!is.na(primary_diagnosis_name) & grepl("cerebrovascular|stroke|tia", primary_diagnosis_name, ignore.case = TRUE), 1, 0),
    
    # Pulmonary conditions  
    charlson_copd = ifelse(has_copd, 1, 0),
    
    # Metabolic conditions
    charlson_diabetes = ifelse(has_diabetes, 1, 0),
    charlson_diabetes_comp = ifelse(!is.na(primary_diagnosis_name) & grepl("diabetes.*complication|diabetic.*retinopathy|diabetic.*nephropathy", primary_diagnosis_name, ignore.case = TRUE), 1, 0),
    
    # Renal conditions
    charlson_renal = ifelse(!is.na(primary_diagnosis_name) & grepl("chronic kidney|renal failure|chronic renal", primary_diagnosis_name, ignore.case = TRUE), 1, 0),
    
    # Liver conditions
    charlson_liver_mild = ifelse(!is.na(primary_diagnosis_name) & grepl("hepatitis|liver disease", primary_diagnosis_name, ignore.case = TRUE) & 
                                !grepl("severe|cirrhosis|failure", primary_diagnosis_name, ignore.case = TRUE), 1, 0),
    charlson_liver_severe = ifelse(!is.na(primary_diagnosis_name) & grepl("liver.*severe|cirrhosis|liver.*failure", primary_diagnosis_name, ignore.case = TRUE), 1, 0),
    
    # Cancer conditions
    charlson_cancer = ifelse(!is.na(primary_diagnosis_category) & primary_diagnosis_category == "Neoplasms" & 
                           (!is.na(primary_diagnosis_name) & !grepl("metastatic|metastasis", primary_diagnosis_name, ignore.case = TRUE)), 1, 0),
    charlson_cancer_mets = ifelse(!is.na(primary_diagnosis_name) & grepl("metastatic|metastasis", primary_diagnosis_name, ignore.case = TRUE), 1, 0),
    
    # Other conditions
    charlson_dementia = ifelse(!is.na(primary_diagnosis_name) & grepl("dementia|alzheimer", primary_diagnosis_name, ignore.case = TRUE), 1, 0),
    charlson_paralysis = ifelse(!is.na(primary_diagnosis_name) & grepl("paralysis|hemiplegia|paraplegia", primary_diagnosis_name, ignore.case = TRUE), 1, 0),
    charlson_peptic = ifelse(!is.na(primary_diagnosis_name) & grepl("peptic ulcer|gastric ulcer", primary_diagnosis_name, ignore.case = TRUE), 1, 0),
    charlson_rheum = ifelse(!is.na(primary_diagnosis_name) & grepl("rheumatoid|lupus|connective tissue", primary_diagnosis_name, ignore.case = TRUE), 1, 0),
    charlson_aids = ifelse(!is.na(primary_diagnosis_name) & grepl("hiv|aids", primary_diagnosis_name, ignore.case = TRUE), 1, 0)
  ) %>%
  mutate(
    # Calculate weighted Charlson score
    charlson_score = 
      (charlson_mi * 1) +
      (charlson_chf * 1) +
      (charlson_pvd * 1) +
      (charlson_cvd * 1) +
      (charlson_copd * 1) +
      (charlson_diabetes * 1) +
      (charlson_diabetes_comp * 2) +
      (charlson_renal * 2) +
      (charlson_liver_mild * 1) +
      (charlson_liver_severe * 3) +
      (charlson_cancer * 2) +
      (charlson_cancer_mets * 6) +
      (charlson_dementia * 1) +
      (charlson_paralysis * 2) +
      (charlson_peptic * 1) +
      (charlson_rheum * 1) +
      (charlson_aids * 6),
    
    # Create risk categories
    charlson_category = case_when(
      charlson_score == 0 ~ "Low Risk (0)",
      charlson_score %in% 1:2 ~ "Moderate Risk (1-2)",
      charlson_score %in% 3:4 ~ "High Risk (3-4)",
      charlson_score >= 5 ~ "Very High Risk (5+)"
    )
  )

charlson_stats = charlson_components %>%
  summarise(
    mean_charlson = round(mean(charlson_score, na.rm = TRUE), 2),
    median_charlson = median(charlson_score, na.rm = TRUE),
    max_charlson = max(charlson_score, na.rm = TRUE),
    q75_charlson = quantile(charlson_score, 0.75, na.rm = TRUE),
    percent_high_risk = round(mean(charlson_score >= 3, na.rm = TRUE) * 100, 1)
  )

charlson_summary_table = data.frame(
  Metric = c('Mean Charlson Score', 'Median Charlson Score', 'Max Charlson Score', 
             '75th Percentile Charlson Score', '% High Risk (Score >=3)'),
  Value = c(charlson_stats$mean_charlson, charlson_stats$median_charlson,
            charlson_stats$max_charlson, charlson_stats$q75_charlson,
            charlson_stats$percent_high_risk)
)

kable(charlson_summary_table, caption = 'Charlson Comorbidity Index Summary')

# Charlson Readmission Analysis
charlson_readmission = charlson_components %>%
  group_by(charlson_category) %>%
  summarise(
    total_cases = n(),
    readmissions = sum(readmit_30day, na.rm = TRUE),
    readmission_rate = round((readmissions / total_cases) * 100, 2),
    .groups = 'drop'
  ) %>%
  arrange(desc(readmission_rate))

kable(charlson_readmission, caption = 'Readmission Rates by Charlson Comorbidity Category',
      col.names = c('Charlson Category', 'Total Cases', 'Readmissions', 'Readmission Rate (%)'))

# Update analysis_data with Charlson features
analysis_data = charlson_components 

# Create healthcare utilization intensity features
utilization_features = analysis_data %>%
  mutate(
    # Length of stay categories
    los_category = case_when(
      los_days < 3 ~ 'Short Stay (<3 days)',
      los_days >= 3 & los_days <= 7 ~ 'Medium Stay (3-7 days)',
      los_days > 7 ~ 'Long Stay (>7 days)',
    ),
    
    # Clinical complexity indicators
    total_clinical_events = total_diagnoses + medication_count + procedure_count + total_lab_tests,
    
    # Polypharmacy risk categories
    polypharmacy_risk = case_when(
      medication_count == 0 ~ 'No Medications',
      medication_count <= 5 ~ 'Low Polypharmacy (1-5)',
      medication_count <= 10 ~ 'Moderate Polypharmacy (6-10)',
      medication_count <= 20 ~ 'High Polypharmacy (11-20)',
      medication_count > 20 ~ 'Very High Polypharmacy (>20)'
    ),
    
    # Laboratory testing intensity
    lab_testing_intensity = case_when(
      total_lab_tests == 0 ~ 'No Labs',
      total_lab_tests <= 10 ~ 'Low Testing (1-10)',
      total_lab_tests <= 50 ~ 'Moderate Testing (11-50)',
      total_lab_tests <= 100 ~ 'High Testing (51-100)',
      total_lab_tests > 100 ~ 'Intensive Testing (>100)'
    ),
    
    # Procedure complexity
    procedure_complexity = case_when(
      procedure_count == 0 ~ 'No Procedures',
      procedure_count <= 2 ~ 'Simple (1-2)',
      procedure_count <= 5 ~ 'Moderate (3-5)',
      procedure_count > 5 ~ 'Complex (>5)'
    ),
    
    # Age-adjusted risk
    age_risk_category = case_when(
      age_at_adm < 30 ~ 'Young Adult(<30)',
      age_at_adm < 50 ~ 'Middle Age (30-49)',
      age_at_adm < 65 ~ 'Older Adult (50-64)',
      age_at_adm < 80 ~ 'Elderly (65-79)',
      age_at_adm >= 80 ~ 'Very Elderly (80+)'
    ),
    
    # High-risk combinations
    high_risk_combination = case_when(
      age_at_adm >= 65 & charlson_score >= 3 ~ 'Elderly + High Comorbidity',
      age_at_adm >= 65 & medication_count > 10 ~ 'Elderly + Polypharmacy',
      charlson_score >= 3 & medication_count > 10 ~ 'High Comorbidity + Polypharmacy',
      age_at_adm >= 80 ~ 'Very Elderly Only',
      charlson_score >= 5 ~ 'Very High Comorbidity Only',
      medication_count > 20 ~ 'Very High Polypharmacy Only',
      TRUE ~ 'Standard Risk'
    )
  )

# Utilization feature statistics
utilization_stats = utilization_features %>%
  summarise(
    mean_clinical_events = round(mean(total_clinical_events, na.rm = TRUE), 1),
    median_clinical_events = median(total_clinical_events, na.rm = TRUE),
    high_polypharmacy_percent = round(mean(medication_count > 10, na.rm = TRUE) * 100, 1),
    intensive_testing_percent = round(mean(total_lab_tests > 100, na.rm = TRUE) * 100, 1),
    high_risk_combo_percent = round(mean(high_risk_combination != 'Standard Risk', na.rm = TRUE) * 100, 1)
  )

utilization_summary_table = data.frame(
  Metric = c('Mean Total Clinical Events', 'Median Total Clinical Events', 
             '% with High Polypharmacy (>10 meds)', '% with Intensive Lab Testing (>100 tests)', 
             '% with High-Risk Combinations'),
  Value = c(utilization_stats$mean_clinical_events, utilization_stats$median_clinical_events,
            utilization_stats$high_polypharmacy_percent, utilization_stats$intensive_testing_percent,
            utilization_stats$high_risk_combo_percent)
)

kable(utilization_summary_table, caption = 'Healthcare Utilization Feature Summary')

# High-risk combination analysis
risk_combo_analysis = utilization_features %>%
  group_by(high_risk_combination) %>%
  summarise(
    total_cases = n(),
    readmissions = sum(readmit_30day, na.rm = TRUE),
    readmission_rate = round((readmissions / total_cases) * 100, 2),
    .groups = 'drop'
  ) %>%
  arrange(desc(readmission_rate))

kable(risk_combo_analysis, caption = 'Readmission Rates by High-Risk Combinations',
      col.names = c('Risk Combination', 'Total Cases', 'Readmissions', 'Readmission Rate (%)'))